*===============================================================================
*Project: Characterizing health provider knowledge: evidence from HIV services 
*		  in Kenya, Rwanda, South Africa, and Zambia

*Objective: Estimate the GLM models described in the paper and reproduce the figures
*Owner: ORPHEA team

*Database required: "final_db_orphea_models.dta"
*===============================================================================

set more off
clear

global dataverse "$DROPBOX\ORPHEA_Quality_2019\Manuscript\Submited paper\Dataverse" //Folder were dataverse data was downloaded


*===============================================================================

use "$dataverse\final_db_orphea_models_dataverse.dta", clear //cambiar nombre

*===============================================================================

*Models

**Naive model**
glm score_Theta i.country##ib2.codes_01_11 i.country##i.intervention i01_01 i01_02 log_px_week years_experience ///
score_total_mgmt i.ownership_2 i.fac_type_2 i.urbanicity imp_urbanicity anual_scale ///
imp_anual_scale age_program_ave imp_age_program_ave dedication_ave imp_dedication_ave i01_09 ///
i01_12 day_person dum_senior_pca, family(poisson) link(identity)
margins, dydx(*) post
est sto Overall_pca

**heteroskedasticity-robust standard errors**
glm score_Theta i.country##ib2.codes_01_11 i.country##i.intervention i01_01 i01_02   log_px_week years_experience ///
score_total_mgmt i.ownership_2 i.fac_type_2 i.urbanicity imp_urbanicity anual_scale ///
imp_anual_scale age_program_ave imp_age_program_ave dedication_ave imp_dedication_ave i01_09 ///
i01_12 day_person dum_senior_pca, family(poisson) link(identity) vce(robust)
margins, dydx(*) post
est sto Overall_pca_robust

** cluster-robust standard errors (facility)**
glm score_Theta i.country##ib2.codes_01_11 i.country##i.intervention i01_01 i01_02   log_px_week years_experience ///
score_total_mgmt i.ownership_2 i.fac_type_2 i.urbanicity imp_urbanicity anual_scale ///
imp_anual_scale age_program_ave imp_age_program_ave dedication_ave imp_dedication_ave i01_09 ///
i01_12 day_person dum_senior_pca, family(poisson) link(identity) vce(cluster ID_facility)
margins, dydx(*) post
est sto Overall_pca_cluster

outreg2 [Overall_pca Overall_pca_robust Overall_pca_cluster] ///
using "$dataverse/interaction_models_last_version.xls", excel replace bdec(3) tdec(3) label 

*===============================================================================

*Fig 3

*(Marginsplot)
glm score_Theta i.country##ib2.codes_01_11 i.country##i.intervention i01_01 i01_02 log_px_week years_experience ///
score_total_mgmt i.ownership_2 i.fac_type_2 i.urbanicity imp_urbanicity anual_scale ///
imp_anual_scale age_program_ave imp_age_program_ave dedication_ave imp_dedication_ave i01_09 ///
i01_12 day_person dum_senior_pca, family(poisson) link(identity)

set scheme sj

	*years in the position*
	margins, at((p10) i01_12) at((p25) i01_12) at((p50) i01_12) ///
	at((p75) i01_12) at((p90) i01_12)
	marginsplot, xtitle("") ytitle("Knowledge score") ylabel(30(1)42, labsize(vsmall) angle(0)) ///
	xlabel(1 "P10" 2 "P25" 3 "P50" 4 "P75" 5 "P90", labsize(small) angle(45)) graphregion(color(white) ///
	icolor(white)) fxsize(36) title("Provider years in position", ///
	size(small)) saving("$dataverse\years_pos_margins.gph", replace)
	
	*Management Score*
	margins, at((p10) score_total_mgmt) at((p25) score_total_mgmt) at((p50) score_total_mgmt) ///
	at((p75) score_total_mgmt) at((p90) score_total_mgmt)
	marginsplot, xtitle("") ytitle("") ylabel(30(1)42, labsize(vsmall) angle(0)) ///
	xlabel(1 "P10" 2 "P25" 3 "P50" 4 "P75" 5 "P90", labsize(small) angle(45)) fxsize(28) ///
	graphregion(color(white) icolor(white)) title("Management score", ///
	size(small)) saving("$dataverse\management_margins.gph", replace)
	
	*Average age program*
	margins, at((p10) age_program_ave) at((p25) age_program_ave) at((p50) age_program_ave) ///
	at((p75) age_program_ave) at((p90) age_program_ave)
	marginsplot, xtitle("") ytitle("") ylabel(30(1)42, labsize(vsmall) angle(0)) ///
	xlabel(1 "P10" 2 "P25" 3 "P50" 4 "P75" 5 "P90", labsize(small) angle(45)) graphregion(color(white) ///
	icolor(white)) fxsize(28) title("Program age", ///
	size(small)) saving("$dataverse\age_program_margins.gph", replace)
	
	*Average proportion of exclusive staff*
	margins, at((p10) dedication_ave) at((p25) dedication_ave) at((p50) dedication_ave) ///
	at((p75) dedication_ave) at((p90) dedication_ave)
	marginsplot, xtitle("") ytitle("") ylabel(30(1)42, labsize(vsmall) angle(0)) ///
	xlabel(1 "P10" 2 "P25" 3 "P50" 4 "P75" 5 "P90", labsize(small) angle(45)) graphregion(color(white) ///
	icolor(white)) fxsize(28) title("Proportion of exclusive staff", ///
	size(small)) saving("$dataverse\exclusive_staff_margins.gph", replace)

graph combine "$dataverse\years_pos_margins.gph" "$dataverse\age_program_margins.gph" "$dataverse\exclusive_staff_margins.gph"  ///
"$dataverse\management_margins.gph", graphregion(color(white) icolor(white)) ycommon rows(1) ///
note("Predictive margins with 95% CIs", size(vsmall)) imargin(0.4 0.4 0.4 0.4) iscale(0.8) saving("$dataverse\fig_3.gph", replace)
graph export "$dataverse\fig_3.tif", width(1800) replace

*===============================================================================

*Fig 1 & 2

clear

use "$dataverse\final_db_orphea_fig1&2.dta", clear

**Figure 1**
twoway (kdensity score_Theta_vig if intervention==1,  range(0 100) lcolor(black) lpattern(shortdash_dot)) ///
(kdensity score_Theta_vig if intervention==2, range(0 100) lcolor(gs5) lpattern(dash)) ///
(kdensity score_Theta_vig if intervention==3, range(0 100) lcolor(gs10)), by(country, cols(1) note("") legend(at(5) pos(0))) ///
ytitle("Density", size(small)) xtitle("Knowledge score", size(small)) ylabel( , labsize(small)) ///
xlabel( , labsize(small)) legend(rows(1) label(1 "HTC") label(2 "PMTCT") label(3 "VMMC") ///
region(lcolor(white))) graphr(ic(white) fc(white) lc(white)) subtitle(, fcolor(gs14)) ///
saving("$dataverse\fig_1.gph", replace)
graph export "$dataverse\fig_1.tif", width(1800) replace

**Figura 2**
twoway (kdensity score_Theta_vig if country==1, range(0 100) lcolor(black) ///
lpattern(shortdash)) (kdensity score_Theta_vig if country==2, range(0 100) ///
lcolor(gs4) lpattern(shortdash_dot)) (kdensity score_Theta_vig if country==3, ///
range(0 100) lcolor(gs8) lpattern(dash)) (kdensity score_Theta_vig if country==4, ///
range(0 100) lcolor(gs12) lpattern(solid)), by(intervention, cols(1) note("") legend(at(4) pos(0))) ///
ytitle("Density", size(small)) xtitle("Knowledge score", size(small)) ylabel( , labsize(small)) ///
xlabel( , labsize(small)) legend(rows(1) label(1 "Kenya") label(2 "Rwanda") label(3 "South Africa") ///
label(4 "Zambia") region(lcolor(white))) graphr(ic(white) fc(white) lc(white)) subtitle(, fcolor(gs14)) ///
saving("$dataverse\fig2.gph", replace)
graph export "$dataverse\fig2.tif", width(1800) replace


